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Abstract. We give accurate estimates for the bond percolation critical probabilities on 
seven Archimedean lattices, for which the critical probabilities are unknown, using an 
algorithm of Newman and Ziff . 

o 

^ 1. Introduction 

< 

Since the introduction of percolation theory, it has been an interesting and challenging 

problem to determine the percolation thresholds. Only a few non-trivial graphs are exactly 

solved, such as the bond model on the square, triangular and hexagonal lattices, and the 

site model on the triangular and Kagome lattices. Recently Scullard and Ziff have found 

many new thresholds for various classes of lattices [SI E]. They also conjecture, [10], the 

value of the threshold for one lattice considered here: the (3, 12 2 ) lattice, and also for the 

Kagome lattice. Precise estimates have been calculated for example for the site model on 

the Archimedean lattices, |7j, and for the bond model on the Kagome lattice, [TT] . 

Bounds, more or less tight, have been found by various authors for some lattices. For the 

bond model on the Archimedean lattices, see [I] for a review of rigorous bounds on critical 

probabilities. Recently, Riordan and Walters, [5], have given tight rigorous confidence 

intervals for both site and bond percolation on all Archimedean lattices. 

In this paper, we provide precise estimates for the bond percolation thresholds for the 

O unsolved Archimedean lattices. 

i O i 

1.1. Archimedean lattices. The Archimedean lattices are the vertex transitive graphs 
that can be embedded in the plane such that every face is a regular polygon. A polygon is 
regular if all edges have the same length, and all interior angles are the same. Kepler, [TJ, 
Q\ showed that there exists exactly 11 such graphs. 

The lattices are given names according to the sizes (number of sides of the polygon) of 
faces incident to a given vertex. The face sizes are listed in order, starting with a face such 
O that the list is the smallest possible in lexicographical order. The square lattice thus gets 

the name (4,4,4,4), abbreviated to (4 4 ), and the Kagome lattice the name (3,6,3,6). 
Square representations of the Archimedean lattices studied here, are shown in Figure [TJ 
The bond percolation threshold is known exactly for the square, triangular, and hexag- 
^ onal lattices (the values are 0.5, 2sin(j|), and 1 — 2sin(y|)), and the threshold for the 

Kagome lattice has previously been estimated to 0.5244053 by Ziff and Suding, [TT]. Ziff 
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and Scullard's, [TU] conjectured value for (3, 12 2 ) is 0.740421178..., and it is consistent 
with the estimate given here: 0.740 4219 (8) — the difference is less than one standard 
deviation (further, there is a small positive bias in the estimate). 

In this work we estimate the thresholds for all the Archimedean lattices with unknown 
(and in one case conjectured) threshold except the previously studied Kagome lattice. The 
hexagonal lattice is used benchmark. 




(4,6,12) (3,12^) (4,8 2 ) 



Figure 1. Finite subgraphs of square embeddings of 8 Archimedean lattices. 

1.2. Percolation. The bond percolation process is defined as follows. For each edge of a 
graph G, declare the edge to be open with probability p, independently of all other edges, 
and closed otherwise. 

For the kind of graphs studied here, it is well known that there exists a critical value 
p c (G), called the critical probability, or threshold, such that for p > p c {G), there exists an 
unique infinite connected component of open edges, while for p < p c {G), only finitely large 
connected components of open edges exist. 

1.3. Motivation. One objective of this simulation study is to get an empirical answer to 
the following question. 

Question 1.1. IfG and H are two Archimedean lattices, with self-avoiding walk connective 
constants fi(G) andfi(H), site percolation thresholds p s c (G) andp s c (H), and bond percolation 
thresholds p b c (G) and p b (H) , is it true that 

KG) < im(H) & p s c (G) > p s c (H) p b c (G) > p b c (H). 
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The answer is negative for general two dimensional, quasi-transitive, planar graphs, as 
shown by Wierman, |S], and Parviainen, [3]. 

Available estimates and rigorous bounds for site percolation thresholds and connective 
constants suggest that these two models give the same order on the Archimedean lattices. 

The estimates given in this work suggest, however, that there exist two pairs of Archimedean 



lattices for which the bond and site percolation thresholds are in opposite order, namely 
the pair (3,6,3,6) (Kagome) and (3,4,6,4) (Ruby), and the pair (3 3 ,4 2 ) and (3 2 ,4,3,4). 



Newman and Ziff, [2], have developed a fast algorithm for estimating percolation thresh- 
olds (both bond and site); the running time is nearly linear in the system size (the number 
of vertices or edges of the subgraph), while still producing accurate estimates. 

It turns out that it is profitable to modify the subgraphs used, and use torus shaped 
regions. We will consider Q(p), the probability that a cluster wraps around the torus in 
one direction, but not both. 

If we generate a percolation process on a subgraph with M edges, by adding edges in 
random order, we get M + l different, but dependent, realizations of a percolation process. 
These give a rough estimate of the function q(n), the probability that a cluster spans the 
torus in one direction but not the other, given that there are exactly n open edges. 

Assume for the moment that q(n) is exactly known. The law of total probability then 
gives us Q(p): 



Thus from a single realisation of the Newman-Ziff algorithm we get an estimate of Q(p) 
for all values of p, from to 1. 

In our case p c can be estimated by the value of p at which Q(p) is maximised, since 
the probability Q(p) tends to zero both for p both above and below p c , as the system size 
grows. 

For the method to achieve its impressive running time, it is necessary to keep track of 
clusters efficiently. This can be achieved by a tree based union/find algorithm. For each 
added edge, we find the clusters to which the endpoints belong. If the clusters are different, 
the union of the clusters is calculated. Both steps are rapidly done by representing the 
set of clusters as a directed forest. By a small modification of the union/find algorithm, 
detection of cluster wrapping is also easy. 

We can speed up the execution by delaying the convolution (equation (|2])) and maximi- 
sation step, by averaging over a batch of, say m, estimates of q(n) and use the average to 
estimate Q(p) and p c . One minor drawback by doing so is that we get fewer samples for es- 
timation the statistical error. Also, the standard erroiQ does not decrease as 1/ sjm — but 
this is only true for small m, and a 1/y/m factor dominates from m w 10 an onwards. For 
the largest system sizes considered here, the standard error for large m is approximately 

1 the standard error is an estimate of the standard deviation of the error of the estimate — if the estimator 
is a mean of n values, the standard error is s.e. — s/y/n where s is the sample standard deviation. 



2. The method of Newman and Ziff 




M-n 
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2/y/m times that for m — 1. We thus lose only a factor 2 which is cheap considering 
that the convolution and maximisation step can be of an order of 100 times more time 
consuming than the generation step. 

We have also run short simulations with the hull gradient method (see for example 
[7]), using different representations (of the covering graphs) of the lattices, to verify the 
implementations. (A reason for not using the hull gradient for longer runs was the memory 
requirement; the simulations were run on standard desktop computers, simultaneously used 
in daily work by the author's colleagues.) 

3. Numerical results 

The finite size error decreases very fast. For our benchmark case, the hexagonal lattice, 
the finite size bias was of order 10~ 8 for systems of size 50 000. 

This prompted us to concentrate the computations on only one large system per lattice. 
We used square shaped subgraphs, with between 60 000 and 75 000 edges, for which we 
believe the finite size errors to be of order 10~ 8 , an order of magnitude smaller than the 
statistical errors. The results are summarized in Table [TJ The standard error given is 
s m /\^n where s m is the sample standard deviation when m iterations are used for each 
maximisation step. If we do n realisations, s m is given by 

iV 1 8=1 

where pW is the estimate from realisation i and p c is the grand estimate, the mean of the 
Approximate (1 — a)% confidence intervals are given by 

p c ± z a / 2 s.e., 

where z a / 2 is the a/2 quantile of the Gaussian distribution^] For all lattices we used batches 
of 10 5 realisations for each maximisation step. 

In [2], it was observed that the standard deviation scales as 

where M is the system size (the number of edges). In our simulations we generally observe 
slightly slower convergence. Fitting a = CM~ a we found some variations in the estimated 
values of a between lattices; we observed values in the range 0.25 to 0.40. The values 
decrease as the number of iterations m used per maximisation step increases, but appear 
to settle down for m ~ 20. The between-lattice variation also decreases with m. As the 
purpose here is to get precise estimates for a number of lattices, we have not studied the 
convergence rate thoroughly for all combinations of lattices and parameters. In Table 1 
we report estimates of a for m = 100. For the hexagonal lattice, Figure 2 shows a log- 
log-log plot of the sample standard deviation as a function of system size and number of 
iterations per maximisation step. The numbers in parenthesis give the estimated values of 
the exponents a for fixed values of m. 

2 Strictly speaking i-distribution quantiles should be used, but for the values of n used here they agree 
with the Gaussian quantiles. 
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Figure 2. Sample standard deviation as function of system size and number 
of iterations per maximisation step. 

Table 1. Simulation results. The standard error is denoted by s.e., the 
number of realizations by n, the system size by M, and the estimated value 
of the exponent a is denoted by a. 



Lattice 


Pc 


s.e./lO" 7 


n 


M 


a 


(3,12") 


0.740 42195 


8.0 


2485 


67500 


0.276 


(4,6,12) 


0.693 733 83 


7.2 


2930 


73 728 


0.277 


(4,8 2 ) 


0.676 802 32 


6.3 


2983 


60 000 


0.284 


(3,4,6,4) 
(3 4 ,6) 
(3 3 ,4 2 ) 
(3 2 ,4,3,4) 


0.524 832 58 


5.3 


7568 


64800 


0.286 


0.434 306 21 


5.0 


7397 


64000 


0.274 


0.419 64191 


4.3 


9822 


64000 


0.274 


0.41413743 


4.6 


7504 


64000 


0.269 
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For the hexagonal lattice we did long simulations for 6 different system sizes to study 
the precision. We used batches of 10 4 realisations. The results are summarized in Table [2j 
The exact value of the bond percolation threshold is p c = 1 — 2 sin(7r/18) ~ 0.652 703 6446. 

Table 2. Simulation results for the hexagonal lattice. Here M denotes the 
number of edges, n the total number of realizations used for the estimate p c , 
with standard error s.e. Further, si is the sample standard deviation when 
one repetition is used for each maximisation step, and is an estimate of 
the limiting sample standard deviation per repetition when m repetitions are 
used for each maximisation step. 



M 


n/10 4 


Pc 


s.e. 




Si 




96 


754 700 


0.652 685 19 


6.68 x 10" 


-7 


0.048 


0.058 


384 


634 300 


0.652 70781 


5.38 x 10" 


-7 


0.032 


0.043 


1536 


460 520 


0.652 704 77 


4.49 x 10- 


-7 


0.021 


0.030 


6144 


444 786 


0.652 703 77 


3.12 x 10- 


-7 


0.013 


0.021 


24 576 


41727 


0.652 703 87 


6.58 x 10- 


-7 


0.0079 


0.014 


98 304 


5326 


0.652 703 67 


1.32 x 10- 


-6 


0.0048 


0.0097 



The small finite size bias observed in [2] is confirmed by our simulations. The estimates 
pc are conjectured to converge like 

\Pc — Pc\ ~ M -11 / 8 . 

Fitting 

V C = p c + bM- ll ' s 

to data, excluding the M = 96 data point, gives b = 0.0017 and a finite size error for 
M = 50 000 of 3.4 x KT 8 . 

Excluding further data points for small values of M gives slight variations, but the finite 
size error is always below 4 x 10 8 . This supports our belief that the finite size error for our 
main estimates are considerable less than the statistical error. 

It is also worthwhile to note that, except for very small sizes, the finite size bias is 
positive, in contrast to the site percolation case studied in [2]. (The same is observed with 
the hull gradient method; for the bond cases studied here, the finite size bias is positive, 
while the site cases on the same lattices studied in [7] have negative finite size biases.) 
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